// -*- C++ -*-

#ifndef CLUSTER_H
#define CLUSTER_H

#include <cmath>
#include <limits>

// The picture channel number

class Point;
class Cluster;

const int CH = 15;

// Point

class Point {
  friend double distance_avg(const Cluster& p, const Point& q);
  friend double mon_distance_avg(const Cluster& p, const Point& q);
public:
  Point() {
    for (int i = 0; i < CH; ++i) {
      data[i] = 0.0;
    }
  }

  bool zero() const {
    for (int i = 0; i < CH; ++i) {
      if (data[i] != 0.0) return false;
    }
    return true;
  }

  const double& operator[](int index) const {
    return data[index];
  }

  double& operator[](int index) {
    return data[index];
  }

  Point& operator*(double p) {
    for (int k = 0; k < CH; ++k) {
      data[k] *= p;
    }
    return *this;
  }

  Point& operator/(double p) {
    for (int k = 0; k < CH; ++k) {
      data[k] /= p;
    }
    return *this;
  }

private:
  double data[CH];
};

inline bool operator<(const Point& p, const Point& q) {
  for (int k = 0; k < CH; ++k) {
    if (p[k] < q[k]) return true;
    if (q[k] < p[k]) return false;
  }
  return false;
}

inline bool operator==(const Point& p, const Point& q) {
  for (int k = 0; k < CH; ++k) {
    if (p[k] != q[k]) return false;
  }
  return true;
}

inline bool operator!=(const Point& p, const Point& q) {
  return !(p == q);
}

// Point - Point distances

double distance_rp(const Point& p, const Point& q, double d = 2.0);
double distance_euclid(const Point& p, const Point& q);
double distance_rinf(const Point& p, const Point& q);
double distance_div(const Point& p, const Point& q);
double distance_jm(const Point& p, const Point& q);

double mon_distance_rp(const Point& p, const Point& q, double d = 2.0);
double mon_distance_euclid(const Point& p, const Point& q);

// Cluster

class Cluster {
  friend double distance_avg(const Cluster& p, const Cluster& q);
  friend double distance_avg(const Cluster& p, const Point& q);
  friend double mon_distance_avg(const Cluster& p, const Point& q);
public:
  
  Cluster() {
    clear();
  }

  void clear() {
    _size = 0.0;
    for (int i = 0; i < CH; ++i) {
      _sum[i] = 0.0; 
      for (int j = 0; j < CH; ++j) {
	_sqr[i][j] = 0.0;
      }
    }
  }

  double size() const {
    return _size;
  }

  double expected(int i) const {
    return _sum[i] / _size;
  }

  double covariance(int i, int j) const {
    return (_sqr[i][j] - _sum[i] * _sum[j] / _size) / _size;
  }

  double deviance(int i) const {
    return std::sqrt(covariance(i, i));
  }

  double avgDeviance() const {
    double r = 0.0;
    for (int k = 0; k < CH; ++k) {
      r += covariance(k, k);
    }
    return std::sqrt(r);
  }

  void scale(double rate) {
    for (int i = 0; i < CH; ++i) {
      for (int j = 0; j < CH; ++j) {
        _sqr[i][j] = (_size * _sqr[i][j] * rate + 
                      _sum[i] * _sum[j] * (1.0 - rate)) / _size;
      }
    }
  }

  Cluster& operator=(const Point& point) {
    _size = 1;
    for (int i = 0; i < CH; ++i) {
      _sum[i] = point[i];
      for (int j = 0; j < CH; ++j) {
	_sqr[i][j] = point[i] * point[j];
      }
    }
    return *this;
  }

  Cluster& operator+=(const Point& point) {
    ++_size;
    for (int i = 0; i < CH; ++i) {
      _sum[i] += point[i];
      for (int j = 0; j < CH; ++j) {
	_sqr[i][j] += point[i] * point[j];
      }
    }
    return *this;
  }

  Cluster& operator+=(const Cluster& cluster) {
    _size += cluster._size;
    for (int i = 0; i < CH; ++i) {
      _sum[i] += cluster._sum[i];
      for (int j = 0; j < CH; ++j) {
	_sqr[i][j] += cluster._sqr[i][j];
      }
    }
    return *this;
  }

  Cluster& operator-=(const Point& point) {
    --_size;
    for (int i = 0; i < CH; ++i) {
      _sum[i] -= point[i];
      for (int j = 0; j < CH; ++j) {
	_sqr[i][j] -= point[i] * point[j];
      }
    }
    return *this;
  }

  Cluster& operator-=(const Cluster& cluster) {
    _size -= cluster._size;
    for (int i = 0; i < CH; ++i) {
      _sum[i] -= cluster._sum[i];
      for (int j = 0; j < CH; ++j) {
	_sqr[i][j] -= cluster._sqr[i][j];
      }
    }
    return *this;
  }

  Cluster& add(const Point& point, double count = 1) {
    _size += count;
    for (int i = 0; i < CH; ++i) {
      _sum[i] += point[i] * count;
      for (int j = 0; j < CH; ++j) {
	_sqr[i][j] += point[i] * point[j] * count;
      }
    }
    return *this;
  }

  Cluster& remove(const Point& point, int count = 1) {
    _size -= count;
    for (int i = 0; i < CH; ++i) {
      _sum[i] -= point[i] * count;
      for (int j = 0; j < CH; ++j) {
	_sqr[i][j] -= point[i] * point[j] * count;
      }
    }
    return *this;
  }

  Point center() const {
    Point center;
    for (int i = 0; i < CH; ++i) {
      center[i] = expected(i);
    }
    return center;
  }

private:
  double _size;
  double _sum[CH];
  double _sqr[CH][CH];
};

inline Cluster operator+(const Cluster& p, const Point& q) {
  Cluster r = p;
  return r += q;
}

inline Cluster operator-(const Cluster& p, const Point& q) {
  Cluster r = p;
  return r -= q;
}

inline Cluster operator+(const Cluster& p, const Cluster& q) {
  Cluster r = p;
  return r += q;
}

inline Cluster operator-(const Cluster& p, const Cluster& q) {
  Cluster r = p;
  return r -= q;
}

// Cluster - Cluster distances

// Quadratic average of the euclidean distance of the points
double distance_avg(const Cluster& p, const Cluster& q);
// Divergence of the clusters
double distance_div(const Cluster& p, const Cluster& q);
// Statistical distance of the clusters
double distance_stat(const Cluster& p, const Cluster& q);
// Mahalanobis Distance
double distance_mah(const Cluster& p, const Cluster& q);
// Bhattacharya Distance
double distance_bat(const Cluster& p, const Cluster& q);

// Cluster - Point distances

// R_p distance
double distance_rp(const Cluster& p, const Point& q, double d = 2.0);
// Average distance
double distance_avg(const Cluster& p, const Point& q);
// Mononton average distance
double mon_distance_avg(const Cluster& p, const Point& q);

// Inverse monoton mappings

// Euclidian distance
double inv_map_rp(double dist, double d = 2.0);
// R_p distance
double inv_map_euclid(double dist);

#endif
